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1 Introduction 

Computer simulations allow for the investigation of many materials properties 
and processes that are not easily accessible in the laboratory. This is particularly 
true in the Earth sciences, where the relevant pressures and temperatures may 
be so extreme that no experimental techniques can operate at those conditions. 
Computer modeling is often the only source of information on the properties 
of materials that, combined with indirect evidence (such as e.g. seismic data), 
allows one to discriminate among competing planetary models. Many computer 
simulations are performed using effective inter-atomic potentials taylored to re- 
produce some experimentally observed properties of the materials being investi- 
gated. The remoteness of the physically interesting conditions from those achiev- 
able in the laboratory, as well as the huge variety of different atomic coordination 
and local chemical state occurring in the Earth interior, make the dependability 
of semi-empirical potentials questionable. First-principles techniques based on 
density-functional theory (DFT) (Hohenberg & Kohn 1964, Kohn & Sham 1965) 
are much more predictive, not being biased by any prior experimental input, 
and have demonstrated a considerable accuracy in a wide class of materials and 
variety of external conditions. The importance of thermal effects in the range of 
phenomena interesting to the Earth sciences makes a proper account of atomic 
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motion essential. Traditionally, this is achieved using molecular dynamics tech- 
niques which have been successfully combined with DFT in the first-principles 
molecular dynamics technique of Car and Parrinello (Car & Parrinello 1985). 
Well below the melting temperature, the numerical efficiency of molecular dy- 
namics is limited by the lack of ergodicity, which would require long simulation 
times, and by the importance of long- wavelength collective motions (phonons) , 
which would require large simulation cells. Both difficulties are successfully 
dealt with in the quasi-harmonic approximation (QHA) where the thermal prop- 
erties of solid materials are traced back to those of a system of non-interacting 
phonons (whose frequencies are however allowed to depend on volume or on 
other thermodynamic constraints) . An additional advantage of the QHA is that 
it accounts for quantum-mechanical zero-point effects, which would not be ac- 
cessible to molecular dynamics with classical nuclear motion. The availability of 
suitable techniques to calculate the vibrational properties of extended materials 
using a combination of DFT and linear-response techniques (resulting in the so- 
called density-functional perturbation theory, DFPT (Baroni et al. 1987, Baroni 
et al. 2001)) makes it possible to combine the QHA with DFT. The resulting 
simulation methodology has shown to be remarkably accurate in a wide tem- 
perature range, extending up to not very far from the melting line and has been 
applied to a wide variety of systems, including many which are relevant to the 
Earth sciences. This paper gives a short overview of the calculation of thermal 
properties of materials in the framework of the QHA, using DFT. The paper is 
organized as follows: in Sec. 2 we introduce some of the thermal properties of 
interest and describe how they can be calculated in the framework of the QHA; 
in Sec. 3 we describe the DFPT approach to lattice dynamics; in Sec. 4 we 
briefly introduce some of the computer codes that can be used to perform this 
task; in Sec. 5 we review some of the application of the first-principles QHA 
to the study of the thermal properties of materials; finally, Sec. 6 contains our 
conclusions. 



2 Thermal properties and the Quasi-Harmonic 
Approximation 

The low-temperature specific heat of solids is experimentally found to vanish 
as the cube of the temperature, with a cubic coefficient that is system-specific 
(Kittel 1996, Wallace 1998). This is contrary to the predictions of classical 
statistical mechanics, according to which the heat capacity of a system of har- 
monic oscillators does not depend on temperature, nor on its spectrum. One 
of the landmarks of modern solid-state physics, that greatly contributed to the 
establishment of our present quantum-mechanical picture of matter, is the De- 
bye model for the heat capacity of solids. This model naturally explains the 
low-temperature specific heat of solids in terms of the (quantum) statistical me- 
chanics of an ensemble of harmonic oscillators, which can in turn be pictorially 
described as a gas of non-interacting quasi-particles obeying the Bose-Einstein 
statistics (phonons). 

The internal energy of a single harmonic oscillator of angular frequency w, 



2 



in thermal equilibrium at temperatiue T, is: 

, huj huj , , 

{E) = ^ + ^ , (1) 

where fc^ is the Boltzmann constant. By differentiating with respect to tem- 
perature the sum over all the possible values of the phonon momentum in the 
Brillouin zone (BZ) of Eq. ((!]), the constant- volume specific heat of a crystal 
reads: 

^^(^) = -^^fiw(q,^K(q,z.), (2) 

where ^(q, v) is the frequency of the i^-th mode (phonon) at point q in the 

BZ, n'(q,!/) = ^ ( e ~ ' ^'^'^ ^^^^ extended to the first BZ. By 

assuming that there are three degenerate modes at each point of the BZ, each 
one with frequency a;(q, v) = c|q|, c being the sound velocity, and converting the 
sum in Eq. (21) into an integral, the resulting expression for the heat capacity, 
valid in the low-temperature limit, reads: 

where fl is the volume of the crystal unit cell and Od = is the 

so-called Debye temperature. 

In the Born-Oppenheimer approximation (Born & Oppenheimer 1927), the 
vibrational properties of molecules and solids are determined by their electronic 
structure through the dependence of the ground-state energy on the coordinates 
of the atomic nuclei (Martin 2004). At low temperature the amplitudes of atomic 
vibrations are much smaller than inter-atomic distances, and one can assume 
that the dependence of the ground-state energy on the deviation from equilib- 
rium of the atomic positions is quadratic. In this, so called harmonic, approxi- 
mation (HA) energy differences can be calculated from electronic-structure the- 
ory using static response functions (DeCicco & Johnson 1969, Pick et al. 1970) 
or perturbation theory (Baroni et al. 1987, Baroni et al. 2001) (See Sec. [3]). 

In the HA, vibrational frequencies do not depend on interatomic distances, 
so that the vibrational contribution to the crystal internal energy does not de- 
pend on volume. As a consequence, constant-pressure and constant-volume 
specific heats coincide in this approximation, and the equilibrium volume of a 
crystal does not depend on temperature. Other shortcomings of the HA include 
its prediction of an infinite thermal conductivity, infinite phonon lifetimes, and 
the independence of vibrational spectra (as well as related properties: elastic 
constants, sound velocities etc.) on temperature, to name but a few. A proper 
account of anharmonic effects on the static and dynamical properties of ma- 
terials would require the calculation of phonon-phonon interaction coefficients 
for all modes in the BZ. Although the leading terms of such interactions can 
be computed even from first principles (Baroni & Debernardi 1994, Debernardi 
et al. 1995) — and the resulting vibrational linewidths have in fact been evaluated 
in some cases (Debernardi et al. 1995, Lazzeri et al. 2003, Bonini et al. 2007) — 
the extensive sampling of the phonon-phonon interactions over the BZ required 
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for free-energy evaluations remains a daunting task. The simplest generalization 
of the HA, which corrects for most of the above mentioned deficiencies, while 
not requiring any explicit calculation of anharmonic interaction coefficients, is 
the QHA. 

In the QHA, the crystal free energy is assumed to be determined by the 
vibrational spectrum via the standard harmonic expression: 

T) = Uo{X) + i ^ fi^(q, y\X) + fesT^] log (l - c^^^^fe^) , (4) 

qi^ qi/ 

where X indicates any global static constraint upon which vibrational frequen- 
cies may depend (most commonly just the volume V ^ but X may also include 
anisotropic components of the strain tensor, some externally applied fields, the 
internal distortions of the crystal unit cell, or other thermodynamic constraints 
that may be applied to the system), and Uq{X) is the zero-temperature energy 
of the crystal as a function oi X . In the case X = V , differentiation of Eq. Q 
with respect to volume gives the equation of state: 

dF 



P 



(5) 



where 



V duj{q,v) 
^(q'")^- ^q,.) dV 
are the so-called Griineisen mode parameters. In a perfectly harmonic crystal, 
phonon frequencies do not depend on the interatomic distances, hence on vol- 
ume. In such a harmonic crystal Eq. ^ implies that the temperature derivative 
of pressure at fixed volume vanish: {dP/dT)v ~ 0. It follows that the ther- 
mal expansivity, j3 = V^^{dV/dT)p, which is given by the therniodynamical 
relation: 

^ JdP/dT)y 

^ idP/dv)T ^ ' 

= ■^^fi^i(l,i^h{(l,'^)n'iq,v), (9) 

^ qi/ 

where BT = V(dP/dV)T is the crystal bulk modulus, would also vanish for 
perfectly harmonic crystals. Inspired by Eq. let us define Cy(q,j^) = 

fia;(q, z^)n'(q, i')/V as the contribution of the f-th normal mode at the q point 
of the BZ to the total specific heat, and 7 as the weighted average of the various 
Griineisen parameters: 

Eq^7(q,i')Cv(q,i') 

In terms of 7, the thermal expansivity simply reads: 

= (11) 
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The vanishing of the thermal expansivity in the HA would also imply the equal- 
ity of the constant-pressure and constant-volume specific heats. By imposing 
that the total differentials of the entropy as a function of pressure and tempera- 
ture or of volume and temperature coincide, and by using the Maxwell identities, 
one can in fact show that (Wallace 1998): 

= TBtP'^. (13) 

We conclude this brief introduction to the QHA by noticing that the ansatz 
given by Eq. ^ for the crystal free energy in terms of its (volume-dependent) 
vibrational frequencies gives immediate access to all the equilibrium thermal 
properties of the system. Whether this implicit account of anharmonic effects 
through the volume dependence of the vibrational frequency only is sufficient to 
describe the relevant thermal effects, or else an explicit account of the various 
phonon-phonon interactions is in order, instead, is a question that can only be 
settled by extensive numeric experience. 



3 Ab-initio phonons 

3.1 Lattice dynamics from electronic-structure theory 

Several simplified approaches exist that allow to calculate full (harmonic) phonon 
dispersions w(q, v) from semi-empirical force fields or inter-atomic potentials 
(Briiesch 1982, Singh 1982). The accuracy of such semi-empirical models is 
however often limited to the physical conditions (pressure, atomic coordination, 
crystal structure, etc.) at which the inter-atomic potentials are fitted. Re- 
ally predictive calculations, not biased by the experimental information used 
to describe inter-atomic interactions require a proper quantum-mechanical de- 
scription of the chemical bonds that held matter together. This can be achieved 
in the framework of electronic-structure theory (Martin 2004), starting from 
the adiabatic or Born and Oppenheimer (BO) approximation, and using mod- 
ern concepts from DFT (Hohenberg & Kohn 1964, Kohn & Sham 1965) and 
perturbation theory (Baroni et al. 2001). 

Within the BO approximation, the lattice-dynamical properties of a system 
are determined by the eigenvalues E and eigenfunctions $ of the Schrodinger 
equation: 

- E + EBoim)^ Hm) = mm), (w) 

where R/ is the coordinate of the I-th nucleus, Mj its mass, {R} indicates the 
set of all the nuclear coordinates, and Ebo is the ground-state energy of a system 
of interacting electrons moving in the field of fixed nuclei, whose Hamiltonian — 
which acts onto the electronic variables and depends parametrically upon {R} — 
reads: 

^^^oim) = -^Y.^ + Y^ + E (ri) + mm), (15) 
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y 2 

— e being the electron charge, V{R}(r) = — |r_Rj is the electron-nucleus 

interaction, and -EAr({R}) = 'Ylii=^] \Rj'^Rj\ inter-nuclear interaction en- 
ergy. The equilibrium geometry of the system is determined by the condition 
that the forces acting on individual nuclei vanish: 

P,,_Ma).o, (16) 

whereas the vibrational frequencies, w, are determined by the eigenvalues of the 
Hessian of the BO energy, scaled by the nuclear masses: 



dct 



1 d^Esoim) 2 

- id 



^/MiMj dRidRj 



0. (17) 



The calculation of the equilibrium geometry and vibrational properties of a 
system thus amounts to computing the first and second derivatives of its BO 
energy surface. The basic tool to accomplish this goal is the Hellmann-Feynman 
(HF) theorem (Hellmann 1937, Feynman 1939), which leads to the following 
expression for the forces: 

f I ^^lW^w ggjv(R) 

where n{R}(r) is the ground-state electron charge density corresponding to the 
nuclear configuration {R}. The Hessian of the BO energy surface appearing 
in Eq. (|17p is obtained by differentiating the HF forces with respect to nuclear 
coordinates: 



d^EBoim) _ d¥i 



9R/9R,7 9Rj 

an{R}(r) 9V{R}(r) 



I 



(19) 

dr (20) 



^^">^^)^R;aR7''+ d^jd^j ■ 

Eq. ((2T|l states that the calculation of the Hessian of the BO energy surfaces 
requires the calculation of the ground-state electron charge density, njR-j.(r), as 
well as of its linear response to a distortion of the nuclear geometry, 9n{R}(r)/9Rj. 
This fundamental result was first stated in the late sixties by De Cicco and 
Johnson (DeCicco & Johnson 1969) and by Pick, Cohen, and Martin (Pick 
et al. 1970). The Hessian matrix is usually called the matrix of the inter-atomic 
force constants (IFC). For a crystal, we can write: 

C,, (22) 

where (R) is the a-th Cartesian components of the displacement of the s-th 
atom of the crystal unit cell located at lattice site R, and translational invariance 
shows manifestly in the dependence of the IFC matrix on R and R' through 
their difference only. 
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3.2 Density-functional perturbation theory 

We have seen that the electron-density linear response of a system determines 
the matrix of its IFCs, Eq. (|2ip. Let us see now how this response can be 
obtained from DFT. The procedure described in the following is usually referred 
to as density-functional perturbation theory (Baroni et al. 1987, Baroni et al. 
2001). 

In order to simplify the notation and make the argument more general, we 
assume that the external potential acting on the electrons is a differentiable 
function of a set of parameters, A = {Ai} (Ai = R/ in the case of lattice 
dynamics). According to the HF theorem, the first and second derivatives of 
the ground-state energy read: 

- ' ^ -n^(r)dr, (23) 



dXi J dX. 

I "n^ir)dr+ J' dr. (24) 



dXidXj J dXidXj J dXi dXj 

In DFT the electron charge-density distribution, n"^, is given by: 

N/2 

(r) = 2^|V'^(r)p, (25) 



=1 



where N is the number of electrons in the system (double degeneracy with re- 
spect to spin degrees of freedom is assumed), the single-particle orbitals, "(/"^(r), 
satisfy the Kohn-Sham (KS) Schrodinger equation: 

+ ^scF(r)) ^^(r) = 6,^^^(r), (26) 

and the self- consistent potential, Vg(jp^ is given by: 

V^CF = V^ + j 0^^dv' + ^ixc[n^]{v). (27) 

where /ijcc is the so-called exchange-correlation (XC) potential (Kohn & Sham 
1965). The electron-density response, dn\{r)/dXi, appearing in Eq. can be 
evaluated by linearizing Eqs. and (P7)l with respect to wave-function, 

density, and potential variations, respectively. Linearization of Eq. (j25p leads 
to: 

N/2 

n'(r)=4Re^C(r)^;W, (28) 

n=l 

where the prime symbol (as in n') indicates differentiation with respect to one 
of the A's. The super-script A has been omitted in Eq. as well as in any 

subsequent formulas where such an omission does not give rise to ambiguities. 
Since the external potential (both unperturbed and perturbed) is real, KS eigen- 
functions can be chosen to be real, and the sign of complex conjugation, as well 
as the prescription to keep only the real part, can be dropped in Eq. 

The variation of the KS orbitals, -(/'^(r), is obtained by standard first-order 
perturbation theory (Messiah 1962): 

(H^SCF - <Wn) = -iV^CF - (29) 



7 



where Hg^jp = — 5^ -§^1 + Vscf unperturbed KS Hamiltonian, , 

V^icF(r) = ^'(r) + j ^(r, v')n\v')dv' (30) 

is the first-order correction to tlie self-consistent potential, Eq. (P7)l . K(r,r') = 
+ is Hartree-plus-XC kernel, and = (Cl^sCFl^n) is the 

first order variation of the KS eigenvalue, e„. Equations ([^51 - 150)) form a set 
of self-consistent equations for the perturbed system completely analogous to 
the KS equations in the unperturbed case — Eqs. (P5|) . (pS)) . and (P7)l — with 
the KS eigenvalue equation, Eq. being replaced by a linear system, Eq. 

(P^l . The computational cost of the determination of the density response to 
a single perturbation is of the same order as that needed for the calculation of 
the unperturbed ground-state density. 

The above discussion applies to insulators, where there is a finite gap. In 
metals a finite density of states occurs at the Fermi energy, and a change in the 
orbital occupation number may occur upon the application of an infinitesimal 
perturbation. The modifications of DFPT needed to treat the linear response 
of metals are discussed in Refs. (de Gironcoli 1995, Baroni et al. 2001). 



3.3 Interatomic force constants and phonon band inter- 
polation 

The above discussion indicates that the primary physical ingredient of a lattice- 
dynamical calculation is the IFC matrix, Eq. (j2ip . from which vibrational 
frequencies can be obtained by solving the secular problem, Eq. ((T7]). That 
phonon frequencies can be classified according to a well defined value of the 
crystal momentum q follows from the translational invariance of the IFC matrix. 
Because of this, the IFC matrix can be Fourier analyzed to yield the so called 
dynamical matrix, prior to diagonalization: 

(q)=^C,f (R)e-'i-^, (31) 

R 

and the squared vibrational frequencies, w(q, z^)^, are the eigenvalues of the 
in X 3n dynamical matrix: 

n being the number of atoms in the unit cell. The direct computation of the 
IFCs is unwieldy because it requires the calculation of the crystal electronic 
linear response to a localized perturbation (the displacement of a single atom 
or atomic plane), which would in turn break the translational symmetry of 
the system, thus requiring the use of computationally expensive large unit cells 
(Martin 2004, Alfe n.d., Parlinski n.d.). The IFCs are instead more conveniently 
calculated in Fourier space, which gives direct access to the relevant q-dependent 
dynamical matrices (Baroni et al. 2001). Because of translational invariance, 
the linear response to a monochromatic perturbation, i.e. one with a definite 
wave- vector q, is also monochromatic, and all quantities entering the calculation 
can be expressed in terms of lattice-periodic quantities (Baroni et al. 2001). As a 
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consequence, vibrational frequencies can be calculated at any wave- vector in the 
BZ, without using any supercells, with a numerical effort that is independent 
of the phonon wave-length and comparable to that of a single ground-state 
calculation for the unperturbed system. 

The accurate calculation of sums (integrals) of lattice-dynamical properties 
over the BZ (such as those appearing in the QHA formulation of the thermody- 
namics of crystals in Sec. [5]) requires sampling the integrand over a fine grid of 
points. This may be impractical in many cases, and suitable interpolation tech- 
niques are therefore called for. The most accurate, and physically motivated, 
such technique consists in the calculation of real-space IFCs by inverse analyzing 
a limited number of dynamical matrices calculated on a coarse grid. Dynamical 
matrices at any arbitrary point in the BZ can then be inexpensively recon- 
structed by Fourier analysis of the IPC's thus obtained. According to the sam- 
pling theorem by Shannon (Shannon 1949), if the IFCs are strictly short-range, 
a finite number of dynamical matrices, sampled on a correspondingly coarse 
reciprocal-space grid, is sufficient to calculate them exactly by inverse Fourier 
analysis. The IFCs thus obtained can then be used to calculate exactly the 
dynamical matrices at any wave- vector not included in the original reciprocal- 
space grid. In the framework of lattice-dynamical and band-structure calcula- 
tions this procedure is usually referred to as Fourier interpolation. Of course, 
IFCs are never strictly short-range, and Fourier interpolation is in general a 
numerical approximation, subject to so-called aliasing errors, whose magnitude 
and importance have to be checked on a case-by-case basis. 

Let us specialized to the case of a crystal, in which lattice vectors R are 
generated by primitive vectors ai , a.2, aa : R/mn = ^^i + 'nia.2 + na.3 , with l,m,n 
integer numbers. The reciprocal lattice vectors G are generated in an analogous 
way by vectors bi, b2, ba, such that 



Correspondingly we consider a symmetry-adapted uniform grid of q-vectors: 



where q, r are also integers. This grid spans the reciprocal lattice of a supercell 
of the original lattice, generated by primitive vectors iViai, N2Sl2: ^^333. Since 
wave- vectors differing by a reciprocal-lattice vector are equivalent, all values of 
p,q,r differing by a multiple of iVi, N2, N3 respectively, are equivalent. We can 
then restrict our grid to p G [0, iVi - 1], g S [0, A^2 - 1], and r € [0, A^a - 1]. 
The qpqr grid thus contains Ni x N2 x uniformly spaced points and spans 
the parallelepiped generated by bi, b2, h^. It is often convenient to identify 
wave- vectors with integer labels spanning the [""j: ~ l] range, rather than 
[0, A^ — 1]. Negative indeces can be folded to positive values using the periodicity 
of discrete Fourier transforms. 

Once dynamical matrices have been calculated on the cihki gi'id, IFCs are 
easily obtained by (discrete) fast-Fourier transform (FFT) techniques: 



a; • b 




(33) 




(34) 



r" 




N1N2N3 



1 
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JVi JV2 
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where the bi-orthogonahty of the real- and reciprocal-space primitive vectors, 
Eq. is used to get q^,^ • R/„„ = Stt ( ^ ^) . The IPCs thus ob- 

tained can be used to calculate dynamical matrices at wave- vectors not originally 
contained in the reciprocal-space grid. This can be done directly wave-vector 
by wave-vector, or by FFT techniques, by padding a conveniently large table of 
IFCs with zeroes beyond the range of those calculated from Fourier analyzing 
the original coarse reciprocal-space grid. 

4 Computer codes 

In order to implement the QHA from first principles, one needs to compute 
the complete phonon dispersion of a crystal for different values of the crys- 
tal volume. This can be done within DFT by the direct or frozen phonon 
method, or by the linear response method (Baroni et al. 2001, Martin 2004). 
The former does not require the use of specialized software beside that needed 
to perform standard ground-state DFT calculations, but is computationally 
more demanding. Some software tools that help analyze the output of stan- 
dard DFT code to produce real-space IFC's and, from these, reciprocal-space 
dynamical matrices are available (Alfe n.d., Parlinski n.d.). As for the linear- 
response approach, two widely known general-purpose packages exist. Quan- 
tum ESPRESSO (Giannozzi et al. 2009) and ABINIT (Gonze et al. 2002). 
In the following we briefly describe the former, as well as another code, QHA, 
that can be used as a post-processing tool to perform QHA calculations start- 
ing from lattice-dynamical calculations performed with many different methods 
(semi-empirical as well as first-principles, frozen-phonon, as well as DFPT). 

4.1 Quantum ESPRESSO 

Quantum ESPRESSO (opEn Source Package for Research in Electronic Struc- 
ture, Simulation, and Optimization) is an integrated suite of computer codes 
for electronic-structure calculations and materials modeling, based on DFT, 
plane waves, pseudopotentials (norm-conserving and ultrasoft) and all-electron 
Projector- Augmented- Wave potentials (Giannozzi et al. 2009) . It is freely avail- 
able under the terms of the GNU General Public License. Quantum ESPRESSO 
is organized into packages. For the purposes of lattice-dynamical calculations 
and QHA applications, the two most relevant ones are PWscf and PHonon. The 
former produces the self-consistent electronic structure and all related computa- 
tions (forces, stresses, structural optimization, molecular dynamics). The latter 
solves the DFPT equations and calculates dynamical matrices for a single wave- 
vector or for a uniform grid of wave- vectors; Fourier interpolation can be applied 
to the results to produce IFCs up to a pre-determined range in real space. The 
effects of macroscopic electric field are separately dealt with using the known 
exact results valid in the long- wavelength limit (Born & Huang 1954). Both 
the electronic contribution to the dielectric tensor, e^o, and the effective charges 
Z* are calculated by PHonon and taken into account in the calculation of in- 
teratomic force constants. Once these have been calculated, phonon modes at 
any wave-vector can be recalculated in a quick and economical way. Anhar- 
monic force constants can be explicitly calculated using the D3 code contained 
in the PHonon package. The volume dependence of the IFCs needed within the 
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QHA is simply obtained numerically by performing several phonon (harmonic) 
calculations at different volumes of the unit cell. 

4.2 The QHA code 

Once the IFC matrix (or, equivalently, the dynamical matrix over a uniform 
grid in reciprocal space) has been calculated, thermodynamical properties can 
be easily calculated using the QHA code (Isaev n.d.). QHA requires in input 
just a few data: basic information about the system (such as atomic masses, 
lattice type) and a file containing IFCs, stored in an appropriate format. QHA 
then calculates and several quantities such as the total phonon density of states 
(DOS), atom-projected DOS, the isochoric heat capacity, the Debye temper- 
ature, zero-point vibration energy, internal energy, entropy, mean square dis- 
placements for atoms, etc. The DOS is obtained via the tetrahedron method 
(Lehmann & Taut 1972), while integrals over the frequency are calculated using 
the Simpson's "3/8 rule". 

5 Applications 

The first investigations of the thermal properties of materials using ab initio 
phonons and the QHA date back to the early days of DFPT theory, when the 
thermal expansivity of tetrahedrally coordinated semiconductors and insulators 
was first addressed (Fleszar & Gonze 1990, Pavone 1991, Pavone et al. 1993). 
Many other applications have appeared ever since to metals, hydrides, inter- 
metallic compounds, surfaces, and to systems and properties of mineralogical 
and geophysical interest. Brief reviews of these applications can be found in 
Refs. (Baroni et al. 2001, Rickman & LeSar 2002); this section contains a more 
up-to-date review, with a special attention paid to those applications that are 
relevant to the Earth Sciences. 

5.1 Semiconductors and insulators 

One of the most unusual features of tetrahedrally coordinated elemental and bi- 
nary semiconductors is that they display a negative thermal expansion coefhcient 
(TEC) at very low temperature. This finding prompted the first applications 
of the QHA to semiconductors, using first a semi-empirical approach (Biernacki 

6 SchefHer 1989), and first-principles techniques in the following (Fleszar & 
Gonze 1990, Pavone 1991, Pavone et al. 1993, Hamdi et al. 1993, Debernardi 
& Cardona 1996, Gaal-Nagy et al. 1999, Rignanese et al. 1996, Xie, Chen, Tse, 
de Gironcoli & Baroni 1999, Eckman et al. 2000, Mounet & Marzari 2005, Zim- 
mermann et al. 2008). The detailed insight provided by the latter allowed one 
to trace back this behavior to the negative Griineisen parameter in the low- 
est acoustic phonon branch and to its flatness that enhances its weight in the 
vibrational density of states at low frequency. This behavior is not observed 
in diamond at ambient conditions — which in fact does not display any negative 
TEC (Pavone et al. 1993, Xie, Chen, Tse, de Gironcoh & Baroni 1999)— whereas 
at pressures larger than ^ 700 GPa the softening of the acoustic Griineisen pa- 
rameters determines a negative TEC. The TEC of diamond calculated in Ref. 
(Pavone et al. 1993) starts deviating from experimental points at T > 600K 
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which was explained in terms of enhanced anharmonic effects at higher tem- 
perature. However, a recent calculation done with a different XC energy func- 
tional (GGA, rather than LDA) (Mounet & Marzari 2005) displayed a fairly 
good agreement with experiments up to T = 1200K, and with results of Monte- 
Carlo simulations (Herrero & Ramirez 2000) up to T = 3000K. Graphite shows 
negative in-plane TEC over a broad temperature range, up to 600K, and the 
calculated TEC for graphene is negative up to 2000K (Mounet & Marzari 2005). 
This is due to a negative Griineisen parameter of the out-of-plane lattice vibra- 
tions along the TM and TK directions (the so called ZA modes, which plays an 
important role in the thermal properties of layered materials, due to the high 
phonon DOS displayed at low frequency because of a vanishing sound velocity 
(Lifshitz 1952, Zabel 2001)). Such an unusual thermal contraction for carbon 
fullerenes and nanotubes was confirmed by molecular dynamics simulations in 
Ref. (Kwon et al. 2004). The heat capacity of carbon nanotubes was calcu- 
lated in Ref. (Zimmermann et al. 2008). The out-of-plane TEC calculated for 
graphite (Mounet & Marzari 2005) is in poor agreement with experiment. This 
is not unexpected because inter-layer binding is mostly due to dispersion forces 
which are poorly described by the (semi-) local XC functionals currently used 
in DFT calculations. 

One of the early achievements of DFT that greatly contributed to its estab- 
lishment in the condensed-matter and materials-science communities was the 
prediction of the relative stability of different crystal structures as a function of 
the applied pressure (Gaal-Nagy et al. 1999, Eckman et al. 2000, Correa et al. 
2006, Liu et al. 1999, Isaev et al. 2007, Mikhaylushkin et al. 2007, Dubrovinsky 
et al. 2007). Thanks to the QHA, vibrational effects can be easily included in 
the evaluation of the crystal free energy, thus allowing for the exploration of the 
phase diagram of crystalline solids at finite temperature. In Refs. (Gaal-Nagy 
et al. 1999, Eckman et al. 2000), for instance, the P — T phase diagram for Si 
and Ge was studied in correspondence to the diamond — > /3-Sn transition. No- 
ticeable changes in the EOS of ZnSe at finite temperature were shown in (Hamdi 
et al. 1993). The phase boundary between cubic and hexagonal BN has been 
studied in Ref. (Kern et al. 1999) using the QHA with an empirical correction 
to account for the leading (explicit) anharmonic effects. Other applications of 
the QHA in this area include the low-temperature portion ot the P — T phase 
diagram for the diamond — )- BC8 phase transition (Correa et al. 2006) and 
the sequence of Rhombohedral (223K) — > Orthorhombic (378K) — > Tetragonal 
(778K) — > Cubic phase transitions in BaTiOa (Zhang et al. 2006) at ambient 
pressure. 

5.2 Simple metals 

The QHA has been widely used to investigate the thermal properties of BCC 
(Quong & Liu 1997, Liu et al. 1999, Debernardi et al. 2001), FCC (Debernardi 
et al. 2001, Li & Tse 2000, Grabowski et al. 2009, Xie et al. 2000, Narasimhan 
& de Gironcoli 2002, Xie, de Gironcoli, Baroni & Scheffler 1999a, Tsuchiya 
2003, Sun et al. 2008), and HOP (Ismail et al. 2001, Ahhoff et al. 1993) metals. 
These works generally report a good agreement with experiments as concerns 
the calculated lattice volume, bulk modulus, TEC, Griineisen parameter, and 
high-pressure/high-temperature phase diagram. Some discrepancies in the tem- 
perature dependence of Cp and TEC might be connected to the neglect of ex- 
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plicit anharmonic effects at high temperatures, as weh as due to overestimated 
cell volumes when using GGA XC functional. In Ref. (Grabowski et al. 2009) 
it was stressed that implicit quasi-harmonic effects dominate the thermal prop- 
erties, being almost two orders of magnitude larger than explicit anharmonic 
ones, irrespective of the XC functional adopted. 

The QHA has also been an important ingredient in the calculation of the 
melting curve of some metals, such as Al (Vocadlo & Alfe 2002), Si (Alfe 
& Gillan 2002), and Ta (Gulseren & Cohen 2002, Taioh et al. 2007), per- 
formed via thermodynamic integration. The vibrational contribution to the 
low-temperature free energy of the crystal phase was shown to be important 
for lighter elements (such as Al), whereas it is negligible for heavier ones, such 
as Ta. The P — T phase diagram for HCP-BCC Mg has been obtained in 
Ref. (Althoff et al. 1993), where it was shown that a proper account of lattice 
vibrations improves the prediction of the transition pressure at room temper- 
ature. Interestingly, in Ref. (Xie, de Gironcoli, Baroni & SchefSer 1999a) it 
was noticed that in the QHA equation of state (EOS) of Ag there exists a 
critical temperature beyond which no volume would correspond to a vanishing 
pressure — thus signaling a thermodynamic instability of the system — and that 
this temperature is actually rather close to the experimental melting tempera- 
ture of Ag. Narasimhan et al. (Narasimhan & de Gironcoli 2002) studied the 
influence of different (LDA and GGA) functionals on the thermal properties of 
Cu. The contribution of lattice vibrations to the phase stability of Li and Sn 
has been studied in Refs. (Liu et al. 1999, Pavone et al. 1998, Pavone 2001): 
a proper account of vibrational effects considerably improves the predictions of 
the low-temperature structural properties of a light element such as Li, which is 
strongly affected by zero-point vibrations (Liu et al. 1999). The large vibrational 
entropy associated with low-frequency modes stabilizes the BCC structure of Li 
(Liu et al. 1999) and /3-Sn (Pavone et al. 1998, Pavone 2001) just above room 
temperature. 

5.3 Hydrides 

One of the best illustrations of the ability of the QHA to account for the effects 
of lattice vibrations on the relative stability of different crystalline phases is pro- 
vided by iron and palladium hydrides, FeH and PdH. FeH was synthesized by 
different experimental groups (Antonov et al. 1980, Badding et al. 1992, Hirao 
et al. 2004) and its crystalline structure was found to be a double hexagonal 
hexagonal structure (DHCP), contrary to the results of ab initio calculations 
(Elsasser et al. 1998) that, neglecting vibrational effects, would rather predict 
a simple HCP structure. The puzzle remained unsolved until free-energy cal- 
culations for FCC, HCP, and DHCP FeH (Isaev et al. 2007) showed that the 
hydrogen vibrational contribution to the free energy actually favors the DHCP 
structure. This is a consequence of the linear ordering of H atoms in HCP FeH, 
which shifts to higher frequencies the mostly H-like optical band of the system, 
with respect to the FCC and DHCP phases. The corresponding increase in the 
zero-point energy makes the DHCP structure — which is the next most favored, 
neglecting lattice vibrations — the stablest structure at low pressure. The quan- 
tum nature of hydrogen vibrations and its influence on the phase stability of 
hydrides was also clearly demonstrated in (Caputo & Alavi 2003, Hu et al. 2007). 
First-principles pseudopotential calculations for PdH have shown that tetrahe- 
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drally coordinated H (B3-type PdH) is energetically favored with respect to oc- 
tahedrally coordinated H (Bl-type PdH), at variance with experimental findings 
(Rowe et al. 1972, Nelin 1971). The quantum-mechanical behavior of hydrogen 
vibrations dramatically affects on the stability of PdH phases at low tempera- 
ture, favoring the octahedral coordination of hydrogen atoms in PdH (Caputo 
& Alavi 2003). As another example, the QHA does not predict any a — ^ /3 
(monoclinic to orthorhombic) phase transition in Na2BeH4 (Hu et al. 2007), 
contrary to the conclusions that were reached from static total-energy calcula- 
tions. Overall, the structural parameters of most alkaline hydrides calculated 
using the QHA turned out to be substantially improved by a proper account of 
zero-point vibrations, both using LDA and GGA XC functionals (more so in the 
latter case) (Roma et al. 1996, Barrera et al. 2005, Lebegue et al. 2003, Zhang 
et al. 2007). 

5.4 Intermetallics 

The QHA has been also successfully applied to the thermal properties of in- 
termetallics and alloys. For example, the Griineisen parameters, isothermal 
bulk modulus, TEC, and constant-pressure specific heat for AlaLi have been 
calculated in (Li & Tse 2000) . The TEC temperature dependence of the tech- 
nologically important superalloys B2 NiAl and LI2 NiaAl, as well as LI2 IrsNb, 
have been studied in Refs. (Wang et al. 2004, Arroyave et al. 2005, Lozovoi & 
Mishin 2003, Gornostyrev et al. 2007). This is a very significant achievement of 
QHA, as it makes it possible very accurate temperature-dependent calculations 
of the misfit between lattice parameters of low-temperature FCC/BCC alloy 
and high-temperature LI2/B2 phases, which plays a considerable role in the 
shape formation of precipitates. It has been found that zero-point vibrations do 
not affect the type of structural defects in B2 NiAl, nor do they change quali- 
tatively the statistics of thermal defects in B2 NiAl (Lozovoi & Mishin 2003). 
Ozolins et al. (Ozolins et al. 1998) and Persson et al. (Persson et al. 1999) have 
studied the influence of vibrational energies on the phase stability in Cu-Au 
and Re-W alloys, using a combination of the QHA and of the cluster-variation 
method. It turns out that lattice vibrations considerably enhance to the stabil- 
ity of CuAu intermetallic compounds and Cu-Au alloys with respect to phase 
separation (Ozolins et al. 1998), as well as to the relative stability of the ordered 
vs. disordered phases at high temperature (Persson et al. 1999). 

5.5 Surfaces 

Ab initio calculations for surfaces coupled with the QHA have been done for 
the past 10 years. For example, an anomalous surface thermal expansion, 
the so called surface pre-melting, has been studied for a few metallic sur- 
faces, such as Al(OOl) (Hansen et al. 1999), Al(lll)(Narasinihan & Scheffler 
1997), Ag(lll)(Xie, de GironcoU, Baroni & Scheffler 19996, Narasimhan & 
Scheffler 1997, Al-Rawi et al. 2001), Rh(OOl), Rh(llO) (Xie & Scheffler 1998), 
Mg(lOTO) (Ismail et al. 2001), Be(lOTO) (Lazzeri & de Gironcoli 2002) and 
Be(OOOl) (Pohl et al. 1998). Hansen et al. (Hansen et al. 1999) noticed that 
the QHA is fairly accurate up to the Debye temperature, above which explicit 
anharmonic effects, not accounted for in this approximation, become impor- 
tant. While no peculiar effects for the surface inter-layer spacing were found for 
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Al(lll) (Narasimhan & Scheffler 1997), for Ag and Rh surfaces it was found 
that the outermost inter layer distance, (ii2, is reduced at room temperature, 
with respect to its bulk value, whereas it is expanded at high temperatures 
(Narasimhan & Scheffler 1997, Xie, de Gironcoh, Baroni & Scheffler 19996, Al- 
Rawi et al. 2001, Xie & Scheffler 1998). The expansion of di2 in the Ag and 
Rh surfaces, as well as in Be(OOOl) (Pohl et al. 1998), is related to the softening 
of some in-plane vibrational modes with a corresponding enhancement of their 
contribution to the surface free energy. Free energy calculations for Be(lOlO) 
(Lazzeri & de Gironcoh 2002) and Mg(lOTO) (Ismail et al. 2001) successfully 
account for the experimentally observed oscillatory behavior of the interatomic 
distances. The large contraction of di2 in Be(lOlO) was explained in terms 
of a strong anharmonicity in the second layer in comparison with the surface 
layer (see also (Marzari et al. 1999)). For Be(OOOl) no oscillatory behavior in 
inter-layer spacings was observed in (Pohl et al. 1998), but an anomalously large 
surface thermal expansion does occur. 

5.6 Earth Materials 

The extreme temperature and pressure conditions occurring in the Earth interior 
make many geophysically relevant materials properties and processes difficult, 
if not impossible, to observe in the laboratory. Because of this, computer sim- 
ulation is often a premier, if not unique, source of information in the Earth 
sciences. By increasing the pressure, the melting temperature also increases, 
so that the temperature range over which a material behaves as a harmonic 
solid is correspondingly expanded, thus making the QHA a very useful tool to 
investigate materials properties at Earth-science conditions. 

Iron, the fourth most abundant element on Earth and the main constituent of 
the Earth core, plays an outstanding role in human life and civilization. In Refs. 
(Kormann et al. 2008, Sha & Cohen 20066, Sha & Cohen 2006a) the thermody- 
namics and thermoelastic properties of BCC Fe have been treated by means of 
the QHA and finite-temperature DFT. The temperature dependence of the cal- 
culated constant-pressure heat capacity deviates from experiment at room tem- 
perature, but a proper inclusion of magnetic effects dramatically improves the 
agreement up to the Curie temperature (Kormann et al. 2008). The calculated 
Debye temperature and low-temperature isochoric heat capacity Cy are in good 
agreement with available experimental data. The magnitude and temperature 
dependence of the calculated Ci2; C44 elastic constants (Sha & Cohen 2006 a) 
are consistent with experiment (Leese & Lord Jr. 1968, Dever 1972, Isaak & 
Masuda 1995) in the temperature range from OK to 1200K at ambient pressure, 
while Cii is overestimated (Sha & Cohen 2006a), likely because of an under- 
estimated equilibrium volume. The ambient-pressure shear and compressional 
sound velocities are consistent with available ultrasonic measurements. The c/a 
ratio of £-Fe has been studied in (Sha & Cohen 2006 c) up to temperatures of 
6000 K and pressures of 400 GPa by using the QHA, resulting in good agree- 
ment with previous calculations (Gannarelli et al. 2003) and X-Ray diffraction 
experiment (Ma et al. 2004). A combination of experiments and calculations 
performed within the QHA was used to show that the FCC and HCP phases of 
nonmagnetic Fe (Mikhaylushkin et al. 2007) can co-exist at very high temper- 
atures and pressures 6600 K and 400 GPa) , due to quite small free-energy 
differences. 
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Bl-type MgO and CaO, MgSiOs perovskite, the aragonite and calcite phases 
of CaCOa, the various polymorphs of aluminum silicate, Al2Si05, silica, Si02 
and alumina, AI2O3 are very important constituents of the Earth's crust and 
lower mantle. Besides, it is believed that the Earth's D" layer is mostly com- 
posed of post-perovskite MgSiOa, while 7-spinel Mg2Si04 is the dominant min- 
eral for the lower part of Earth's transition zone. Note that Mg-based minerals 
do contain some amount of Fe substituting Mg. The high-pressure crystalline 
structure and stability of these phases are discussed in (Oganov 2004, Oganov, 
Price & Scandolo 2005). Lattice dynamics and related thermal and elastic 
properties of Bl MgO have been studied in Refs. (Strachan et al. 1999, Drum- 
mond & Ackland 2002, Oganov et al. 2003, Oganov & Dorogokupets 2003, Karki 
et al. 2000, Karki et al. 1999, Wu et al. 2008, Wu & Wentzcovitch 2009). Wentz- 
covich and et al. have introduced a semi- empirical ansatz that allows for an 
account of explicit anharmonic contributions to the QHA estimate of various 
quantities, such as the TEC and Cp (Wu et al. 2008, Wu & Wentzcovitch 2009), 
resulting in a much improved agreement with experiments. The temperature 
and pressure dependence of elastic constants of Bl MgO (Karki et al. 2000, Karki 
et al. 1999, Isaak et al. 1990) calculated within QHA show very good agreement 
with experimental data (Isaak et al. 1989). Besides, pressure dependence of ab 
initio compressional and shear sound velocities is in consistent with seismic 
observations for the Earth's lower mantle (Karki et al. 1999). In contrast with 
these successes, the calculated thermal properties of the Bl and B2 phases of 
CaO (Karki & Wentzcovitch 2003) are inconsistent with experimental data, and 
this is most likely due to the too small lattice parameter predicted by the LDA, 
as later investigations based on a GGA XC functional seem to indicate (Zhang 
& lai Kuo 2009). 

The thermal properties of MgSiOs and Mg2Si04 and the phase transition 
boundary in these minerals (perovskite— >-post-perovskite MgSiOa and spinel— ^-post- 
spinel Mg2Si04) have been extensively studied(Wentzcovitch et al. 2006, Oganov 
& Ono 2004, Oganov k Price 2005, Wu & Wentzcovitch 2009, Yu et al. 2007, Wu 
et al. 2008, Yu et al. 2008, Ono & Oganov 2005) due to their great importance 
for the Earth's D" layer and lower mantle, respectively. Improved EOS of Bl 
MgO (Wu et al. 2008), obtained by means of renormalized phonons and QHA, 
has been used as a new pressure calibration to re-evaluate the high pressure 
- high temperature phase boundary in MgSi03 and Mg2Si04 minerals using 
experimental data from (Fei et al. 2004, Hirose et al. 2006, Speziale et al. 2001). 

Alumina, AI2O3, plays an important role in high-pressure experiments: for 
example, it serves as a window material for shock-wave experiments. Cr-doped 
alumina, ruby, is used as a pressure calibration material in diamond-anvil-cell 
experiments. Besides, it is a component of solid solutions with MgSiOa poly- 
morphs that have significantly different thermal properties from pure MgSiOs 
minerals. Corundum (q;-A1203) is the most stable phase of alumina at ambient 
conditions, preceded by the 6 phase at lower temperature. The energy differ- 
ence between the 9 and a phases of alumina is rather small, and this raised a 
question as to whether a-Al203 is stabilized by phonons. Zero-point vibrations 
stabilize the corundum phase at low temperatures (Lodziana & Parlinski 2003), 
whereas free-energy calculations show that the a phase can not be stabilized 
by phonons only at room temperature. QHA calculations revealed that at high 
pressures alumina transforms to Calr03- (Oganov & Ono 2005) and U2S3-type 
(Umemoto & Wentzcovitch 2008) polymorphs. 
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The P — T phase diagram for Al2Si05 polymorphs (andalusite, siUimanite, 
and kyanite) (Winkler et al. 1991) and the thermal properties of CaCOa poly- 
morphs (calcite and aragonite) (Catti et al. 1993, Pavese et al. 1996) have been 
studied within the QHA using model inter-atomic potentials. The effect of zero- 
point vibrations on the equilibrium volume in the calcite phase was found to 
be quite important and actually larger than the thermal expansion at relatively 
high temperature (Catti et al. 1993). These calculations (Pavese et al. 1996) 
were not able to account for the experimentally observed (Rao et al. 1968) neg- 
ative in-plane TEC in calcite. The heat capacity and entropy calculated for 
the aragonite phase substantially deviate from experiment. All these problems 
can be possibly traced back to the poor transferability of model inter-atomic 
potentials. 

The thermal properties of the a-quartz and stishovite phases of Si02 have 
been studied in (Lee & Gonze 1995). The heat capacities of both phases were 
found to be in good agreement with experimental data (Lord & Morrow 1957, 
Holm et al. 1967), with the stishovite phase having a lower capacity below 
480K. Interestingly, zero-point vibration energy of the stishovite phase affects 
on thermodynamical properties stronger than in the a-quartz phase (Lee & 
Gonze 1995). The P ~ T phase diagram of Si02 has been examined in Refs. 
(Oganov, Gillan & Price 2005, Oganov & Price 2005), with emphasis on the 
stishovite— >CaCl2^a-Pb02—>pyrite structural changes, resulting in a sequence 
of transitions that do not correspond to any observed seismic discontinuities 
within the Earth. Further investigations at ultrahigh temperature and pressure 
show that Si02 exhibits a pyrite— ^-cotunnite phase transition at conditions that 
are appropriate for the core of gas giants and terrestrial exoplanets (Umemoto 
et al. 2006). 

6 Conclusions 

The QHA is a powerful conceptual and practical tool that complements molec- 
ular dynamics in the prediction of the thermal properties of materials not too 
close to the melting line. In the specific case of the Earth Sciences, the QHA 
can provide information on the behavior of geophysically relevant materials at 
those geophysically relevant pressure and temperature conditions that are not 
(easily) achieved in the laboratory. Large-scale calculations using the QHA 
for geophysical research will require the deployment of a large number of re- 
peated structure and lattice-dynamical calculations, as well as the analysis of 
the massive data generated therefrom. We believe that this will require the use 
of dedicated infrastructures that combine some of the features of massively par- 
allel machines with those of a distributed network of computing nodes, in the 
spirit of the grid computing paradigm. The Quantum ESPRESSO distribu- 
tion of computer codes is geared for exploitation on massively parallel machines 
up to several thousands of closely coupled processors and is being equipped with 
specific tools to distribute lattice-dynamical calculation over the grid (Di Meo 
et al. n.d.). 
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